library(tidyverse)
library(readr)
library(ProbBayes)

oby_data <- 5
sd_data <- 1

posterior_density <- function(theta){
  
  prior <- dcauchy(theta, location=0, scale = 1)  # Prior density
  
  likelihood <- dnorm(oby_data, mean=theta, sd=sd_data)  # Likelihood density
  
  return(prior * likelihood)  # Posterior density
}

metropolis_cauchy <- function(n=50000, start=0, scale=1, proposal_sd=1){
  
  samples <- numeric(n)
  samples[1] <- start
  
  for(xi in 2:n){
    current <- samples[xi - 1]
    
    # Generate a proposal from a normal distribution
    proposal_sample <- rnorm(1, mean=current, sd=proposal_sd)
    
    # Calculate the acceptance ratio
    ratio <- min(1, posterior_density(proposal_sample) / posterior_density(current))
    
    # print(ratio)
    
    # Accept or reject the proposal
    if(runif(1) < ratio){
      samples[xi] <- proposal_sample
    } else {
      samples[xi] <- current
    }
  }
  
  return(samples)
}

samples <- metropolis_cauchy(n=50000, start=0, scale=1, proposal_sd=1)


# 后验样本的直方图
hist(samples, breaks = 100, probability = TRUE,
     main = "Posterior Samples (Cauchy Prior + Normal Likelihood)",
     xlab = expression(theta), col = "skyblue")

# 可选：叠加核密度
lines(density(samples), col = "red", lwd = 2)

# 可选：轨迹图
plot(samples, type = "l", col = "gray", main = "Trace Plot", ylab = expression(theta))


metropolis_f2 <- function(logpost, current, C, iter, ...){
  S <- rep(0, iter)
  n_accept <- 0
  
  for(j in 1:iter){
    candidate <- runif(1, min= current - C, max = current + C)
    
    prob <- exp(logpost(candidate, ...) - logpost(current, ...))
    
    accept <- ifelse(runif(1) < prob, "yes", "no")
    
    current <- ifelse(accept == "yes", candidate, current)
    
    S[j] <- current
    
    n_accept <- n_accept + ifelse(accept == "yes", 1, 0)
  }
  
  list(S=S, accept_rate = n_accept / iter)
}

lpost <- function(theta, s){
  dcauchy(theta, s$loc, s$scale, log = TRUE) +
    dnorm(s$ybar, mean=theta, sd=s$se, log = TRUE)
}

buffalo <- read_csv("data/buffalo_snowfall.csv")

data <- buffalo[59:78,c('SEASON', 'JAN')]

ybar <- mean(data$JAN)
print(nrow(data))
print(sd(data$JAN))
se <- sd(data$JAN) / sqrt(nrow(data))

s <- list(loc = 10, scale = 2, ybar = ybar, se = se)

post_samples <- metropolis_f2(lpost, current = 5, C = 20, iter = 10000, s = s)






